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Abstract 

In aerospace computational fluid dynamics calculations, high aspect ratio, or stretched, 
triangulations are often used to adequately resolve the features of a viscous flow around 
bodies. In this paper, we explore alternatives to the Delaunay triangulation which can be 
used to generate high aspect ratio triangulations of point sets. The method is based on a 
variation of the lifting map concept which derives Delaunay triangulations from convex 
hull calculations. 


1 Introduction 

In computational fluid dynamics (CFD) applications, the problem domain must be discretized 
into meshes (or grids) over which the governing equations of fluid dynamics are solved. In 
general, calculations which take the viscosity of fluids into account result in very high solution 
gradients normal to surfaces, and very small solution gradients tangent to surfaces, with the 
effect diminishing with distance from the surface. Effective use of the computational effort 
required to solve the governing flow equations results when a grid contains high aspect ratio, 
surface- conforming elements near bodies and low aspect ratio elements at a distance. This 
study investigates more robust alternatives to the primarily heuristic methods currently used 
to generate high aspect ratio triangulations. 

2 Approaches in Use 

Heuristic methods have been primarily used to generate high aspect ratio elements in un- 
structured grids. These approaches either combine structured grid generation methodologies 
to generate points followed by a triangulation phase to create the grid, or generate the points 
and create the triangulation in tandem. Delaunay triangulation is very often unsuitable for 
the first approach; examples of the inadequacies are described in [1, 2]. Alternative trian- 
gulations have been used, such as the min-max triangulation, mainly because they appear 
to generate grids better suited to CFD applications. Locally optimal, rather than globally 
optimal, methods are used to generate such grids, primarily due to the cost involved. An- 
other purely heuristic approach is the advancing front technique [3, 4]. This technique uses 
a greedy method to generate points. 
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Figure 1 : A distance function resulting in an invalid triangulation 


A heuristic method to generate stretched triangulations was used by Mavriplis [ 5 ]. He 
recast the 11-dimensional planar triangulation problem to an (11+I )-dimensional surface tri- 
angulation problem. A control surface was introduced which was related to the aspect ratio, 
or amount of stretch, desired. Points from the plane were projected to the control surface, 
and Delaunay triangulation was performed 011 the surface. Since the method was intimately 
tied to the point generation method, an assumption of local planarity was valid. 

3 A Robust Alternative 

A key concept behind the Delaunay triangulation is the empty circumcircle criteria. The 
method used by Mavriplis essentially modified the concept of distance 011 the plane used 
to define the empty circumcircle. By assuming local planarity, circumcircles on the control 
surface project to circumellipses on the plane. This can be viewed as related to the use of 
a convex distance function, here an ellipse, to generate a triangulation, similar to the idea 
of convex distance functions for Voronoi diagrams as described by Chew and Drysdale [ 6 ]. 
The technique used by Mavriplis remains heuristic, however, because of the local planarity 
assumption. To be truly general, the control surface should be independent of the point 
generation method, and a general surface triangulation should be employed. A sparse dis- 
tribution of points projected to a control surface with rapidly changing surface gradients 
could result in an invalid triangulation; Figure 1 illustrates how a rapidly changing distance 
function in combination with a sparse set of points could produce overlapping triangles. To 
generalize the convex distance function formulation to be applicable to the problem at hand, 
a distance function must vary throughout the plane. Sufficient care must be taken to insure 
that invalid triangulations could not arise. 

To develop an approach which yields high aspect ratio triangulations that is more robust 
than the ad hoc methods currently in use, we turn to the formulation that relates Delaunay 
triangulation to convex hulls. A point set lifted onto a paraboloid via orthogonal projection 
generates a convex hull; when the lower hull is projected back to the plane, this results 
in a Delaunay triangulation [ 7 ]. It is possible to construct other convex bodies in which 
the convex hull of a point set lifted onto the body similarly produces a triangulation; these 
triangulations then show a “bias” corresponding to the predominant circumshapes. 

First consider the equations of a paraboloid and an arbitrary plane which cuts through 
the paraboloid. From 
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In other words, by eliminating z we get an expression in x and y which is the equation 
for a circle, or the projection of the intersection of the plane in 3-space and the paraboloid 
on the z = 0 or 2D plane. 

Now, consider a convex object f(x,y). To determine the circumshapes on the plane 
centered at (xo,yo) which correspond to the intersection of a plane z = ax + by + c and 
the convex object f(x,y ), we pose the problem as follows: For any point xq and yo on the 
2D plane, we find the plane in 3-space which is tangent to the convex body at f(xo,yo). 
The intersection of the tangent plane and convex object is simply the point f(xo,yo), which 
projects to the point (xo,yo). To “see” the circumshapes centered about this point, the 
tangent plane must be “pushed” into the convex body by an amount 8 to get more than the 
trivial intersection. 

So, the circumshapes about a point (xo,yo) are defined by: 
f(x,y) - (ax + by + c) 
a 

b 
c 

The discussion section describes the results of tests run for different functions f(x,y). 
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4 An alternative approach 

Another way in which the triangulation/convex hull approach might be modified is by chang- 
ing the projection method used to lift the points to the body. Perspective projection along 
the line connecting the projection point (0, 0, z pro j), (x,y, 0), and a convex body is used in- 
stead of orthogonal projection. For this section, we use the paraboloid f(x,y ) = x 2 + y 2 , 
although other convex bodies could be used, as in the previous section. With this method, it 
is possible under certain conditions for projected points to “miss” the paraboloid entirely, so 
this scenario is to be avoided by insuring that z pro j is sufficiently distant. The circumshapes 
are defined by: 
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The method used to determine circumshapes on the 2D plane is similar to that used in the 
previous section; here, the actual values for the partial derivatives are shown since f(x,y ) is 
known in this case. A sketch of the derivation of the circumshape function is given in section 

A. 

An alternative approach is to define a new surface which is based on the perspective 
projection. Orthogonal projection to this new surface is used so the original (x,y) points 
are unchanged, but the corresponding z value obtained is what would have resulted from 
true perspective projection from a point (0, 0, z pro j ) through (x, y, 0) to a paraboloid (in this 
case, z pro j < 0). The equation of such a surface follows; a sketch of its derivation is given in 
section B. 
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5 Discussion of Figures 

In this section we discuss several examples of the variation in circumshape for triangulations 
derived from either convex body alternatives to the paraboloid, or alternative projection. 
The circumshapes display the various biases of the functions which may be exploited in the 
triangulations. 

Figure 4 shows a set of 200 points randomly selected from a uniform distribution within 
the box with corners at ( — 5, —5), (5, —5), ( — 5,5), and (5,5). Each subsequent example uses 
this test data. It may have been scaled or shifted for different test cases, but any modification 
to the test set, along with the reason for it, will be indicated during discussion. 

The circumshapes derived from both orthogonal and perspective projection examples are 
shown for regularly spaced (xo,yo) for each f(x,y ), except as noted in Figure 12. These in 
general do not correspond to a convex hull, but do illustrate the predominant circumshapes in 
a particular portion of the 2D plane which result from planes slicing through the convex sur- 
face. Since the functions used are continuous, one can infer smoothly changing circumshapes. 

The triangulation results in Figures 5-12 are obtained by lifting the test set to the test 
body f(x,y ) using orthogonal projection, Ending the lower convex hull with respect to the 
2D plane, and projecting these results back to the 2D plane. The triangulation results in 
Figure 13 are obtained by lifting the test set to a paraboloid using perspective projection, 
finding the lower convex hull with respect to the projection point (0, 0, z pro j ), and projecting 
these results back (via perspective projection) to the 2D plane. 
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Figure 5 illustrates results for a paraboloid, f(x, y) = x 2 + y 2 . This result corresponds to 
Delaunay triangulation. 

Figure 6 illustrates results for f(x , y) = x 2 + 10 y 2 ] circumshapes are plotted for 8 = 0.05. 
A little manipulation of the equations show that the circumshapes about a point ( xq , yo ) are 

|x j x ° |2 +(</-») 2 = < 

which describes ellipses, with centers depending on ( xo,yo ), and identical eccentricity and 
orientation throughout the plane. 

Figure 7 illustrates results for f(x,y ) = x 2 + y 4 ; circumshapes are plotted for 6 = 0.02. 
The shapes are constant along one coordinate direction. 

Figure 8 illustrates results for f(x,y ) = x 4 + y 4 ; circumshapes are plotted for 6 = 0.02. 
Circumcircles predominate along x = y and x = —y, with increasing elongation towards the 
coordinate axes. 

Figure 9 illustrates results for f(x,y ) = x 3 + y 3 ; circumshapes are plotted for 8 = 0.09. 
The predominant circumshape orientation changes over the plane. Because this body is not 
convex over the entire plane, the data was translated to the upper right quadrant where the 
portion of the body to which points are lifted would be convex. 

Figure 10 illustrates results for f(x,y ) = a; 1-5 + y 1 ' 5 ', circumshapes are plotted for 8 = 
0.015. As in the previous example, the data was translated to the upper right quadrant. 

Figure 11 illustrates results for f(x,y ) = a; 1-5 + 10 y 1 ' 5 ', circumshapes are plotted for 
8 = 0.015. The data was translated to the upper right quadrant. Although the shapes 
appear constant along one coordinate direction, compression of the circumshapes occurs 
from row to row. 

Figure 12 illustrates results for f(x,y ) = \ ^ ^ 4 °° 2 — ^—400 + J 2 °°°2 ^ , a surface de- 
rived from perspective projection considerations with z pro j = —100. The circumshapes 
shown are derived from the convex hull corresponding to the triangulation; a sampling of 
circumshapes are shown for clarity. Because this convex body was not defined for all points 
on the 2D plane, the data set was scaled by a factor of | in both coordinate directions. Note 
that the alternative approach results in circumshapes similar in orientation and eccentricity 
to the true perspective approach which follows. 

Figure 13 illustrates results for perspective projection; z pro j = —100, 8 = 0.05. Because 
the projection from z pro j = —100, was not defined for all points, the data set was scaled 
by a factor of | in both coordinate directions. The circumshapes exhibit radial variation in 
orientation, with eccentricity increasing for increasing distance from the origin. 

6 Conclusion 

We have demonstrated how a modification of the lifting map formulation which derives a 
triangulation from the convex hull of a convex body can be used to produce alternatives to 
the Delaunay triangulation. By suitable construction of the convex body, triangulations can 
then be “biased” to produce stretched or high aspect ratio triangulations. The appeal of this 
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method is that a valid triangulation of the point set will always exist for a corresponding 
convex body, and the method can be generalized to higher dimensions. 
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Figure 2: The relationship between points on the plane and on the paraboloid (in cross-section). 

(a) z p > 0; (b) z v < 0 


A Derivation of Perspective Projection Circumshape Func- 
tion 


Given the equations for a paraboloid, z' = x' 2 -\-y' 2 , a plane, z' = ax' + by ! + c, and a projection 
point (0, 0, z p ), the circumshapes on the 2D plane can be determined in the following manner, 
lysing similar triangles from the diagrams in Figure 2, the following relationships hold: 
Figure 2a: 

~p_ _ x _ y 
z' x — x' y — y' 


Figure 2b: 


Both cases yield: 
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Substituting the expressions for x' and y' in the plane equation yields: 


z' = ax l--)+by l--)+c 


or 


2 = z n 


ax + by + c 


ax + by + z p 

for ax + by + z p y 0. Substituting the expressions for x' and y' in the paraboloid equation 
yields: 

z' = x 2 (l--) 2 + y 2 (l--) 2 


z' = (x 2 + y 2 )(l--Y 


rearranging gives 


z' = (x 2 + y 2 )(l — 2— + —^) 

z v z l p 

z'z 2 v = (x 2 + y 2 )(z' 2 - 2 z'zp + z 2 p ) 
z' 2 (x 2 + y 2 ) - z'(2zp(x 2 + y 2 ) + z 2 ) + z 2 (x 2 + y 2 ) = 0 


'■‘p ) 1 ~p\ 

Substitute the equation for z' from the plane equation into the last result: 

( ax + by + c \ 2 2 2 , ax + by + c ( 2 2 2 \ , 2 , 2 , 2^ n 

\z p — — (x + y ) — z p — — ( 2 Zp(x +y ) + z )+z (x +y) = 0 

V ax + by + z,p J ax + by + z, p \ l > 1 

Expanding and rearranging gives: 

(x 2 + y 2 )( 1 ) - (ax + by + c)( h — + 1) = 0 
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Figure 3: The 2D problem of determining x' . 


B Derivation of Body from Perspective Projection 

The goal is to derive an expression for a body f(x,y) where orthogonal projection is used 
with a s value which would have been obtained from perspective projection. 

The symmetry of the paraboloid about the s-axis can be exploited to derive an expression 
which performs as required. 

First consider the 2D perspective projection problem shown in Figure 3, where z pro j < 0. 
The goal is to find an expression f(x) = x' such that g{x') = x r2 , or g(f(x)) = x 12 . From 
consideration of similar triangles: 


— v 


proj 


x' 2 — z. 


proj 


Rearranging in terms of x' gives: 


12 , ~proj / n 

X H ; -X ~ Z„ ro j = 0 
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Solving for x' yields f(x): 


fix) = x' = 
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proj 


and g(f(x)) is 


g(f(x)) = x r = - 
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Because the paraboloid is circularly symmetric about the s-axis, it may be observed that 
by looking at planes through the s-axis, the 3D problem is identical to the 2D problem. 
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The goal is to find an expression h(x,y ) = z' . In the 2D problem, x is the distance from 
the z-axis. For a 3D formulation, first find the distance to the point under consideration on 
the 2D plane, which is \f x 2 + y 2 . Then use this new “ad” in the 2D formulation (which is 
also equivalent to decomposing this value into its x' and y' components and calculating the 
paraboloid value). This substitution yields 


z = 


1 
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\Jx 2 + r 
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proj 
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proj 


x 2 + y 2 


In this case, the negative square root produces the correct behavior for the function, i.e. , 
a convex object. 
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Figure 4: Test data used for all examples. 
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Figure 13a: Circumsliapes predicted from perspective projection, z pro j = —100, 8 = 0.05 



Figure 13b: Triangulation derived from perspective projection, z pro j = —100 
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